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1.  INTRODUCTION 

This  report  describes  three  new  features  which  have  been  Incorporated  In 
the  EPIC-2  [1]  code: 

1.  Subcycling,  which  permits  different  time  steps  to  be  used  in 
different  parts  of  the  mesh. 

2.  A  four-node  quadrilateral  which  uses  one  quadrature  point  and  employs 
a  consistent  control  of  the  hourglass  modes. 

3.  A  rigid  interface  which  can  be  used  to  connect  meshes  with  different 
elements  sizes. 

As  has  been  succinctly  summarized  by  Zukas  [2]  in  his  survey  of  computer 
methods  for  Impact  and  penetration,  "computational  techniques  have  advanced  to 
the  point  where  extremely  difficult  situations  can  be  analyzed  quickly  .  .  . 
one  dimensional  problems  In  tens  of  minutes  to  a  few  hours,  and  three  dimen¬ 
sional  problems  In  a  few  hours  to  tens  of  hours."  With  the  decreasing  costs 
of  computers,  these  figures,  which  apply  to  the  larger  main-frame  computers, 
promise  to  diminish  rapidly.  Nevertheless,  It  Is  quite  clear  that  we  are 
still  far  from  the  stage  where  even  two-dimensional  calculations  can  be  used 
effectively  In  the  design  and  decision-making  process,  since  the  parametric 
studies  which  are  essential  In  such  highly  nonlinear  simulations  simply  cannot 
be  made  within  the  normal  framework  of  engineering  time  schedules  because  of 
their  long  running  times. 

However,  substantial  savings  in  computer  time  can  be  achieved  in  these 
programs  through  improvements  of  the  time  integration  process  and  element 
technology.  Even  relatively  simple  techniques  such  as  those  described  here 
can  lead  to  reductions  of  computer  time  by  factors  of  as  much  as  3  to  10  for 
two-dimensional  calculations,  reductions  as  large  as  10  to  20  for  three- 


dimensional  calculations.  These  savings  promise  to  be  Important  not  only  In 
the  use  of  such  programs  on  large  ma1n*frame  computers,  but  on  smaller 
computers,  where  they  may  make  these  calculations  feasible. 

In  the  next  section,  we  review  alternative  Integration  methods  and  two 
methods  of  multi -time-step  explicit  time  Integration  (explicit-explicit 
partition)  and  describe  our  reasons  for  the  choice  of  the  method  that  has  been 
programmed  Into  the  EPIC-2  code.  We  then  describe  the  actual  programming 
steps  In  more  detail  and  the  procedures  for  using  this  multi -time  step 
logic.  The  explicit-explicit  algorithm  permits  an  almost  arbitrary  number  of 
different  time  steps  to  be  used  throughout  the  mesh,  as  long  as  the  ratio  of 
time  steps  between  adjacent  elements  are  Integer  multiples  of  each  other.  It 
also  permits  the  time  steps  to  be  varied  and  calculated  automatically,  which 
Is  a  development  which  we  have  not  seen  reported  for  any  other  program. 

A  4-node  quadrilateral  element  [3]  has  been  added  In  EPIC-2  to  provide  an 
alternative  to  the  triangular  element.  Four-node  quadrilaterals  with  one- 
point  quadrature  provide  more  efficient  computations  because  they  require 
fewer  constitutive  evaluations  and  converge  faster  because  they  are  not  as 
stiff  as  triangles;  triangular  elements  tend  to  lock  for  Incompressible 
materials  unless  special  arrangements  of  elements  are  used.  However,  the  one- 
point  quadrature  element  possesses  spurious  singular  modes  known  as  hourglass 
modes.  To  control  these  modes,  a  consistent  hourglass  control  procedure  first 
described  In  Ref.  [3]  Is  used.  Both  this  element  and  the  existing  triangular 
element  are  programmed  to  accommodate  the  explicit-explicit  time  Integration 
scheme.  The  Implementation  of  this  quadrilateral  element  Is  stated  In  Section 
4. 

The  rigid  Interface  provides  a  technique  for  combining  meshes  with 
different  element  sizes  without  any  Intermediate  elements.  It  enhances  the 


efficiency  of  the  subcycling  technique  because  it  increases  the  differences  in 
the  critical  time  steps  of  different  element  groups.  When  previous  versions 
of  the  EPIC-2  mesh  generator  were  used,  the  critical  time  step  of  different 
element  size  groups  did  not  vary  much,  so  subcycling  was  quite  ineffective. 
This  rigid  interface  is  described  in  Section  5. 

2.  RATIONALE  FOR  CHOICE  OF  TIME  INTEGRATION  METHOD 

We  will  write  the  governing  equations  for  the  finite  element 
semi discretization  in  the  form 

M  u  =  f  (2.1) 

where  M  is  an  assembled  mass  matrix  which  is  considered  to  be  lumped  and  hence 
diagonal,  u  is  the  matrix  of  nodal  displacements,  superposed  dots  denote  time 
derivatives,  and  f  is  a  column  matrix  of  nodal  forces  which  is  given  by 


(2.2) 


In  the  central  difference  method  with  a  changing  time  step,  these  equations 
are  integrated  by  updating  the  nodal  velocities  and  displacements  with  the 
following  formulas 


(2.5) 


aF*  »  y  (At”  +  At”"^) 


where  At  Is  the  time  Increment  and  superscripts  denote  the  step  number.  The 
time  increment  At  is  limited  for  stability  by  the  requirement 
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where  la  is  the  maximum  natural  frequency  in  the  finite  element  mesh 
max 

and  uj  is  the  fraction  of  critical  damping  in  this  frequency.  The  formulas  of 
Flanagan  and  Belytschko  [3]  may  be  used  to  show  that  the  maximum  frequency  of 
a  uniform  strain  element  is  bounded  by 
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where  x  and  y  are  the  Lame  constants,  p  the  density,  N  is  the  number  of  nodes 
for  each  element,  A  the  area  of  the  element  and  j  are  the  components  of  the 
strain-displacement  matrix  which  relate  the  velocity  gradients  to  the  nodal 
velocities  by 
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Although  in  [3]  this  formula  is  only  given  in  this  form  for  quadrilaterals,  it 
has  been  recently  shown  that  it  also  applies  to  triangles.  By  using  the  stan¬ 
dard  expressions  for  the  elements  of  the  B  matrix,  and  for  the  dilatational 
wave  speed,  which  is 


4 


(2.9) 
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it  can  be  shown  that  this  formula  gives  the  following  constraint  on  the 


maximum  frequency 


> 

max 
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(2.10) 


where  s  is  a  characteristic  length  of  the  element.  For  a  triangle,  s  is  the 


length  of  the  longest  side.  Using  the  Rayleigh’s  theorem,  it  can  be  shown 


that  the  maximum  mesh  frequency  is  always  bounded  from  above  by  the  maximum 


element  frequency.  For  an  undamped  system,  Eqs.  (6)  and  (10)  can  be  shown  to 


bound  the  time  step  by 


^^crit  ^  sc. 


(2.11) 


which  for  a  right-triangle  gives 


min 

^^crit  ^ 


(2.12) 


where  h^^^  is  the  distance  from  the  node  at  the  right  angle  to  the  hypotenuse. 


Equation  (12)  differs  from  the  formula  used  in  EPIC-2  by  the  factor  2 


in  the  denominator.  This  formula  is  only  a  lower  bound,  so  it  cannot  be  said 


with  certainty  that  a  computation  which  omits  this  factor  would  necessarily  be 


unstable,  but  our  work  indicates  that  in  general  the  omission  of  this  factor 


may  lead  to  any  unconservative  estimate  of  the  stable  time  step. 


Nevertheless,  the  important  point  which  is  made  by  Eqs.  (11)  and  (12)  is 


that  the  stable  time  step  depends  on  the  dimension  of  the  element  and  de- 


creases  as  the  size  of  the  element  decreases.  Therefore  in  penetration  cal¬ 
culations,  where  the  elements  adjacent  to  the  penetrator  are  severely  com¬ 
pressed,  the  stable  time  step  is  drastically  reduced.  This  leads  to  the  main 
drawback  of  single-step,  explicit  integration  where  the  entire  mesh  must  then 
be  integrated  with  this  very  small  time  step. 

Several  new  methods  have  been  devised  for  circumventing  this  difficulty 
of  single  time  step  integration.  Belytschko  and  Mullen  [4],  [5]  have  presen¬ 
ted  an  explicit-implicit  partition  where  part  of  the  mesh  is  integrated 
implicitly,  and  the  rest  is  integrated  explicitly.  They  have  shown  that  the 
stability  limit  on  the  time  step  is  then  determined  strictly  by  the  highest 
frequency  in  the  explicit  partition,  or  in  other  words,  that  Eq.  (6),  and  the 
resulting  Eqs.  (11)  and  (12),  need  only  be  met  in  the  part  of  the  mesh  that  is 
integrated  explicitly. 

At  first  glance,  this  would  seem  to  offer  substantial  benefits  in  the 
penetration  problem,  for  by  integrating  the  highly  compressed  zones  next  to 
the  penetrator  implicitly,  the  severe  reduction  in  the  stable  time  step  which 
is  caused  by  the  squashing  of  these  elements  would  be  avoided.  However,  the 
difficulty  in  applying  this  to  the  penetration  problem  is  that  the  zones  which 
are  compressed  are  adjacent  to  the  slide  line.  Integrating  these  nodes 
implicitly  has  the  following  drawbacks: 

1)  it  is  difficult  to  formulate  a  stiffness  matrix  for  the  slide  line 
because  of  the  constant  realignment  of  nodes  which  takes  place; 

2)  slide  lines  tend  to  malfunction  if  a  large  time  step  is  used. 

An  alternative  explicit-implicit  partition  has  been  developed  by  Hughes 
and  Liu  [6]  and  recently  extended  by  Liu  and  Belytschko  [7]  to  multi -time  step 
explicit-implicit  partitions,  where  different  time  steps  can  be  used  in  con¬ 
junction  with  the  mesh  partitioning.  This  would  appear  to  have  substantial 


potential  in  the  penetration  problem  In  that  the  Implicit  integrator  could  be 
used  in  all  zones  close  to  the  penetrator  except  for  those  directly  on  the 
slide  line.  The  elements  connected  to  the  slide  line  could  then  be  integrated 
with  a  very  small  time  step  and  an  explicit  integrator.  Although  this  alter¬ 
native  is  quite  appealing,  it  was  ruled  out  for  the  following  reasons: 

1)  we  have  had  very  little  experience  with  the  application  of  this 
method  in  structural  dynamics;  the  problems  which  have  been  tried  by  this 
method  have  been  primarily  heat  conduction  problems  which  tend  to  be 
inherently  more  stable; 

2)  the  Implementation  of  this  procedure  would  require  the  development  of 
a  stiffness  matrix  and  hence  considerable  recoding. 

The  extensive  recoding  which  would  be  required  makes  this  alternative 
quite  unattractive.  The  code  structure  of  explicit-implicit  codes  is  inher¬ 
ently  quite  different  from  purely  explicit  codes,  so  that  the  introduction  of 
the  implicit  option  without  complete  recoding  would  be  very  difficult. 

Another  alternative  to  overcoming  the  drawbacks  of  conditional  stability 
is  to  use  different  time  steps  in  different  parts  of  the  mesh.  This  procedure 
was  studied  in  [8],  where  it  was  shown  that  if  linear  interpolation  is  used  on 
the  displacements  of  the  interface  nodes,  the  procedure  is  conditionally 
stable  provided  Eq.  (6)  is  satisfied  for  each  node  and  is  the  lowest 
maximum  frequency  of  any  element  connected  to  the  node.  In  [9]  it  was  shown 
that  this  procedure  is  equivalent  to  a  constant  velocity  interpolation,  which 
is  easier  to  program;  [10]  describes  the  implementation  of  this  method  in  the 
code  SAMSON  2. 

Explicit/explicit  partitions  have  also  been  discussed  by  Wright  [11]  but 
no  details  on  the  implementation  were  given.  For  these  reasons,  the  explicit/ 
explicit  partition  first  investigated  in  [8]  and  further  developed  in  [9]  and 
[10]  was  implemented  in  EPIC-2. 


-I* 


.|*»^»»? .44' , 


3.  EXPLICIT-EXPLICIT  PARTITION  IN  EPIC-2 


The  procedure  used  in  modifying  EPIC-2  is  based  on  an  explicit-explicit 
partitioning  procedure,  or  subcycling  procedure,  presented  in  [10].  In  this 
method,  the  elements  are  separated  into  element  groups,  each  of  which  can  be 
integrated  with  a  different  time  step  subject  to  the  following  restrictions: 

1.  All  time  steps  must  be  integer  multiples  of  the  smallest  time  step. 

2.  If  any  node  is  shared  by  elements  in  two  different  integration 
groups,  the  time  steps  in  these  groups  must  be  integer  multiples  of 
each  other. 

The  time  step  for  each  element  group  is  recomputed  at  the  end  of  the 
total  cycle  but  kept  constant  during  the  subcycles.  All  elements  near  the 
slave  nodes  must  be  included  in  element  group  1  to  ensure  that  the  slide  line 
is  always  integrated  by  the  minimum  time  Increment.  If  the  smallest  time  step 
does  not  belong  to  element  group  1,  the  run  stops  automatically.  No 
additional  requirement  for  the  arrangement  of  element  group  numbers  is 
necessary.  The  elements  adjacent  to  the  master  nodes  can  be  marked  by  group 


number  other  than  1.  As  soon  as  the  slide  line  interaction  is  detected,  the 
group  indicator  for  the  elements  involved  will  be  checked  and  designated  to  be 
1.  This  is  equivalent  to  an  expanding  interaction  zone  in  the  target  and  thus 
avoids  integrating  the  elements  that  are  not  engaged  in  interaction  by  an 
unnecessarily  small  time  step. 

For  purpose  of  defining  how  the  explicit-explicit  partitioning  procedure 
works,  we  will  define  the  following  variables. 


8 


NTGRP:  number  of  groups  into  which  the  finite  element  mesh  is  subdivided 

Atg:  the  time  increment  for  element  group  G 

t^n  :  the  master  time 

At^:  the  master  time  increment,  which  corresponds  to  the  minimum  Atg 
among  all  element  groups  G 

As  stated  previously,  all  element  time  step  increments  Atg  must  be 

integer  multiples  of  At^.  The  maximum  At^  is  called  At 

M  G  max 

The  program  is  designed  so  that  it  automatically  decides  the  appropriate 
nodal  time  step.  This  is  accomplished  by  using  the  largest  time  increment  for 
any  element  group  connected  to  the  node  to  update  the  node.  In  order  to  pro¬ 
gram  this  algorithm,  each  node  therefore  requires  two  additional  words  of 
storage:  the  nodal  time  t^^  and  the  time  increment  used  for  that  node,  Atj^^  . 

The  essence  of  the  procedure  is  as  follows.  We  call  the  time  steps 
necessary  to  advance  the  master  clock  by  At^^^^  a  cycle.  Within  a  cycle,  when¬ 
ever  the  master  clock  t|^|  is  incremented  by  At^^,  all  elements  are  first 
checked.  If  any  element  is  in  a  group  which  is  not  ahead  of  the  master  time, 
i,e.  if  element  I  is  in  group  G  and 

tg  <  t^  (3.1) 

that  element  is  updated.  This  updating  involves  the  calculation  of  new 
velocity  strains,  stresses  and  internal  forces.  The  element  internal  forces 
are  then  added  into  the  global  internal  force  matrix. 

After  all  elements  have  been  updated,  the  nodal  loop  for  updating 
velocities  and  displacements  is  executed.  In  this  loop,  prior  to  updating  the 
velocities  and  displacements,  the  nodal  clock  t^^  is  checked.  If 


i  *  '.a-A*  a-AL 


<  tn  (3.2) 

the  nodal  clock  is  behind  the  master  clock,  so  the  node  is  updated.  In 
addition,  the  nodal  clock  is  updated  using  the  time  increment  for  that  node. 

The  algorithm  assumes  that  a  velocity  strain  formulation  is  used  for  all 
element  calculations.  When  an  element  needs  to  be  updated,  the  latest  avail¬ 
able  velocity  is  used  to  compute  the  velocity  strain.  This  means  if  an  ele¬ 
ment  is  connected  to  a  node  with  a  larger  time  step,  it  uses  the  same  nodal 
velocity  for  all  intermediate  time  steps.  This  corresponds  to  a  constant  ve¬ 
locity  interpolation  or  a  linear  displacement  interpolation,  which  experience 
has  shown  to  be  stable.  A  flow  chart  for  the  procedure  is  shown  in  Table  1. 

An  important  characteristic  of  impact/penetration  calculations  is  that 
tne  nodes  on  the  sliding  interface  cannot  be  updated  with  a  Courant  number  of 
1.  Instead  a  substantially  smaller  time  step,  such  as  0.2  to 

0.5  should  be  used.  This  restriction  is  necessary  because  the  velocity 

adjustment  algorithms  which  are  used  to  enforce  compatibility  on  a  slideline 
fail  if  the  nodes  penetrate  too  far  into  an  element  during  a  time  step. 
Therefore,  in  EPIC-2  the  elements  with  any  nodes  on  the  slide  line  are 
automatically  allocated  to  the  element  group  with  the  smallest  time  step. 
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Table  1 


Flowchart  of  Explicit-Explicit  Partition  { 

1.  Set  initial  condition  ! 

I 

I 

I 

u°  =  u(o)  ,  y."  “  u(o)  j 

i 

! 

initial  accelerations  are  assumed  to  vanish  »  0.  ! 

•<« 

2.  Initialize  clocks  and  cycle  counters 

master  time 

for  all  element  groups  G 
for  all  nodes  N 


3.  Set  up  nodal  time  step  and  subcycle  counter 

n2  =0 

At|^  =  max(AtQ)  for  all  nodes  N  .  G  represents  any  element  group 
G 

connected  to  node  N. 

4.  Update  nodes  behind  master  clock 

a.  00  N  »  1  to  NNODE 

b.  if  t|^  >  t|^,  skip  node  N 

c.  new  accelerations  u^  =  f 

d.  update  nodal  velocities:  u =  u^"  +  atj^  u’jjj 

e.  update  nodal  displacements:  u^^^  “  !iN 

f.  update  nodal  clock: 


t„  =  0  ; 

^G  =  0 

t^  =0 

n  =0 


Compute  Internal  forces 

d  •  Z0PO  I 

b.  DO  N  =  1  to  NELE 


c.  if  tg  >  t^,  skip  element  N  (element  N  belongs  to  group  G) 

d.  compute  velocity  strains:  &  u 

e.  compute  stress  increments  the  constitutive  relations 

f.  update  stress:  ‘‘‘ 


g.  compute  element  internal  forces:  N  ~  ^ 

h.  if  n2  =  1,  compute  stable  time  increment  for  element 

i.  assemble  f*)^}  „  into  f?^} 

~int,N  ~lnt 

j.  update  element  group  clock:  fg  fg  +  ^fg 


compute  f' 


update  master  clock  and  cycle  counters 
tn  t^  +  At„ 
n  n  +  1 


n2  n2  +  1 


8.  if  n2  =  n2  ,  set  new  element  group  time  increment  At^,  n2  ,  =1  and  go 

max  3  r  g»  iiiax 

to  3;  otherwise  go  to  4. 


4.  CONSTANT  STRAIN  QUADRILATERAL  WITH  HOURGLASS  CONTROL 


An  underintegrated  4-node  quadrilateral  element  with  orthogonal  hourglass 
control,  proposed  in  [3J  by  Flanagan  and  Belytschko,  has  been  added  to  the 
EPIC-2  program.  The  element  is  adapted  for  both  plane  strain  and  axisymmetric 
cases.  It  not  only  reduces  the  number  of  elements  by  half  against  a 
triangular  element  mesh  with  the  same  amount  of  nodes,  but  also  permits  a 


greater  stable  time  step 


The  related  equations  for  this  element  can  be  found  in  Ref.  [3].  A 
concise  and  efficient  computational  procedure  Is  summarized  by  Belytschko  et 
a1.  in  [12]. 

In  the  new  version  of  EPIC-2,  the  implementation  of  the  constant  strain 
quadrilateral  element  is  as  following: 

1.  Set  up  the  discrete  gradient  operator  B  such  that  j  =  Nj  ^  (evaluated 
at  the  quadrature  point).  The  lower  case  subscript  runs  from  1  to  2  and 
represents  r  and  z  spatial  coordinates,  respectively.  The  upper  case 
subscript  has  a  range  of  1  to  4  and  represents  the  nodes  of  an  element. 


’^24 

^31 

^42 

^13 

-''42 

'*13 

'*24 

''31' 

(4.1a) 


^  IJ  “  ^  I  ■  ^  J 


''  IJ  "  I  ■  '’d 


^  “  ^’’31  ^42  "  ^'42  ^31^^^ 


(4.1b) 

(4.1c) 

(4. Id) 


2.  Form  the  velocity  gradient 


“i.J  “  "l.J  “11  ■  <“ll  ■  “l3>®jl  *  <“12  -  “l4>®j2 
Note  that  =  -Bjj  and  B^^  =  -Bj2* 


(4.2) 
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3.  Designate  the  mean  radius  for  an  element 


r  »  (Aj  Pj  +  Ag  r2)/(Aj  +  A2) 


where 


\  "  ^21  ^31  "  ^31  ^21 


^  “  ^31  ^41  “  ^31 


?!  »  (r^  +  Tg  +  r3)/3 


^2  “  ^'*1  ^  '*3  ■*■  ^4^^^ 


4.  Calculate  the  strain  and  spin  rate 


*  1  /  *  •  \ 


U..=-i(u..-U.  .)  1#j 

ij  2  '  j,1  1,j^  '' 


and  the  additional  strain  component 


(4.3a) 

(4.3b) 

(4.3c) 

(4.3d) 

(4.3e) 

(4.4a) 

(4.4b) 


■ee 


=  (I 


iipl)/(4r) 


for  axisymmetric  case  only.  (4.5) 


5.  Compute  and  based  on  constitutive  law  and  Von  Mises 

yield  criterion.  This  step  remains  unchanged  as  It  Is  in  subroutine 
STRESS  of  the  original  EPIC-2  program 


where 


h^  »  [1  -1  1  -1]  (4.6b) 

and  the  hourglass  strain  rate 

'’I  <"•') 

7.  Update  the  hourglass  stress 

«l  ■  “  ('®jj  “jj  “ll  -  “14  '5'“’ 

here  k  is  a  user-controlled  parameter  which  determines  how  much  resistance 
would  be  added.  The  recommended  range  for  <  is  from  0.01  to  0.1. 

8.  Compute  the  nodal  force 

(4.9) 

for  plane  strain  element  V  =  A  and  a..  *  0; 
for  ax i symmetric  element  V  =  ZitFA 


Program  Implementation 

Two  additional  subroutines,  QSTRN  and  QFINT,  were  programmed  In  the 
latest  version  of  EPIC-2,  They  are  corresponding  to  the  calculation  of 
velocity  strains  and  Internal  forces  for  the  quadrilateral  element 
respectively. 

The  major  modifications  for  the  explicit-explicit  subcycling  Integration 
scheme  were  made  In  subroutine  LOOP.  We  will  first  define  the  primary 
variables  In  these  modifications,  then  define  what  size  these  arrays  need  to 
be  In  Table  2. 

Major  Variables 


DTN00( N)  . 

.  .time  increment  for  node  N,  Atj^ 

OTNOOO(N)  . 

.  .  previous  time  Increment  for  node  N 

CLKNOO(N)  . 

.  .  clock  time  for  node  N,  t^^ 

DTGRP(NG)  . 

p 

.  .time  Increment  for  element  group  N  , 

^^NG 

DTGRPO(NG). 

.  .  previous  time  Increment  for  element 

group  NG 

CLKGRP(N)  . 

.  .  clock  time  for  element  group  N** 

NEGRP( N)  . 

.  .  Integration  group  number  for  element 

N;  If  NEGRP{N) 

element  N  will  be  Integrated  with  time  Increment  At 

TIME 

master  time  tn^j 

NCYC2  = 

subcycle  counter,  n2  In  flow  chart 

N2LIM  = 

number  of  subcycles  In  cycle 

NTGRP  = 

number  of  element  groups 

NNODE*  = 

number  of  nodes 

NELE*  = 

number  of  elements 

*  already  exist  In  program 


Table  2 


Minimum  Array  Size  In  COMMON  Statement 
of  Variables 


Variable 

Minimum  Size 

DTNOl) 

NNOOE 

OTNODO 

NNOOE 

CLKNOD 

NNOOE 

DTGRP 

NTGRP 

OTGRPO 

NT6RP 

NEGRP 

NELE 

CLKGRP 

NTGRP 

5.  RIGID  INTERFACE  ALGORITHM 

In  order  to  conserve  computational  effort,  the  mesh  generator  in  EPIC-2 
has  the  capability  of  increasing  the  element  size  for  those  elements  away  from 
the  domain  of  Interest,  namely,  the  impact  zone.  This  mesh  gradation  as  it  is 
Implemented  In  EPIC-2  Is  Illustrated  In  Fig,  la;  region  A  represents  the  do¬ 
main  of  Interest.  For  a  uniform  material,  this  type  of  mesh  gradation  results 
In  little  difference  In  the  critical  time  step  among  the  element  groups  since 
the  shortest  side  of  the  element  governs  the  time  step  and  it  is  of  approxi¬ 
mately  the  same  length  In  the  three  groups  of  elements.  The  advantage  of  the 
explicit  subcycling  scheme,  which  depends  on  the  differences  on  critical  time 
steps  between  element  groups,  is  therefore  quite  minimal  with  this  element 


gradation  scheme. 

To  take  advantage  of  subcycling,  it  is  necessary  to  construct  meshes  with 
gradations  in  element  size  that  retain  nearly  square  elements.  Unfortunately, 
this  requires  a  layer  of  transition  elements  between  element  groups  as  shown 
in  Fig,  lb.  However,  this  usually  causes  severe  inconvenience  in  generating 
data.  To  avoid  these  difficulties,  a  rigid  interface  algorithm  was 
developed.  Two  dissimilar  meshes  are  allowed  to  meet  on  the  rigid  interface 
as  shown  in  Fig.  Ic. 

We  designate  the  nodes  on  the  coarse  mesh  side  of  the  rigid  interface  as 
primary  nodes  and  those  on  the  fine  mesh  side  as  secondary  nodes.  The  primary 
nodes  and  secondary  nodes  are  marked  by  P  and  S  respectively  in  Fig.  2. 
Sometimes  a  primary  node  and  a  secondary  node  may  share  the  same  position  in 
the  space,  however  they  are  considered  distinctive  nodes  in  this  algorithm. 

The  ground  rules  of  this  algorithm  are: 

1)  the  primary  nodes  can  move  independently; 

ii)  the  secondary  nodes  must  respond  accordingly,  in  a  manner  which  minimizes 
the  violation  of  compatibility,  to  the  movement  of  the  adjacent  primary 
nodes. 

The  secondary  nodes  are  constrained  to  remain  on  the  line  connecting  two 
adjacent  primary  nodes.  Full  compatibility  is  then  maintained  only  if  the 
ratios  of  the  sides  of  adjacent  elements  are  integer  multiples.  However,  when 
the  ratio  is  not  an  integer  mutiple,  as  in  Fig.  Ic,  the  deviations  from 
compatibility  are  insignificant. 

The  procedure  for  the  rigid  interface  algorithm  is  as  follows: 

Step  1.  Set  up  the  non-dimensional ized  natural  coordinates  for  all 
secondary  nodes.  We  define 
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I 

I 


I  xl  -  xl  I 

~5  ~r 


1  »  1  to  NS 


(5.1) 


where  NS  is  the  number  of  secondary  nodes  between  primary  nodes  and  P^, 

X  represents  the  position  vector  for  a  node  and  i  t  denote  the  magnitude  of  a 
vector.  Subscripts  P  and  S  refer  to  primary  and  secondary  nodes, 
respectively.  Superscripts  denote  the  node  number.  An  example  for  the 
natural  coordinate  is  illustrated  in  Fig.  2. 

Step  2.  Distribute  the  lumped  mass  of  the  secondary  nodes  to  the  adjacent 
primary  nodes. 


NS  . 

*  mi  +  I  (1  -  Cjm’ 

K  K  i^l  1  i 

.  .  NS  . 

Mp  =  <  +  I  C. 


(5.2a) 

(5.2b) 


where  mj,  and  m^  are  the  original  masses  from  the  assembly  operation,  and  Mp 
is  the  adjusted  mass  for  primary  node  P^ . 

Remark :  Steps  (1)  and  (2)  a.'e  executed  before  entering  the  integration  loop 
since,  as  will  be  seen  later,  the  natural  coordinates  remain  constants 
throughout  the  run. 

Step  3.  Compute  the  accelerations,  velocities  and  displacements  of  all  the 
nodes  except  the  secondary  nodes  in  the  regular  manner.  The  adjusted  mass 
Mp,  instead  of  the  original  mass  m’,  must  be  used  for  updating  the  primary 
node  velocity,  i.e. 


(xi) 


(5.3) 
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Step  4.  Enforce  the  rigid  interface  constraint  on  the  velocity  based  on 
1  7 

Ip  and  Vp  of  the  adjacent  primary  nodes  and  ^2’ 
i  1  2 

=  (1  -  C^-)Vp  +  Vp  (loop  over  secondary  nodes  i)  (5.4) 

Then  the  displacement  d^  is  updated  in  the  usual  way. 

Step  5.  No  modification  is  needed  in  the  element  calculations.  After  all 
nodal  forces  are  determined,  the  nodal  forces  on  secondary  nodes  are 
transferred  onto  the  adjacent  primary  nodes. 


f^  +  fj  +  (1  -  5.)fi 

~p  ~p  1  ~S 


(5.5a) 


(loop  over  secondary  nodes  i) 


f^^  ^  f^  +  5. 

~p  ~p  1  ~S 


(5.5b) 


Since  the  displacement  between  two  adjacent  primary  nodes  is  linear  and  a 
linear  interpolation  is  used  in  step  (4)  to  determine  the  secondary  node 
velocity,  the  natural  coordinates  for  secondary  nodes  will  remain  constants. 
This  corresponds  with  our  remark  in  step  (2). 

Two  runs  with  the  rigid  interface  have  been  examined.  The  original  and 
deformed  models  are  shown  in  Figs.  3.  and  4.  No  numerical  instability  has 


been  detected  in  either  case. 
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Figure  3b.  Exam 


6.  INPUT  CHANGES 


The  following  additions  in  input  are  necessary  to  run  this  modified  version  of  EPIC-2. 

Solid  Material  Cards  (6E16.8)  -  Add  one  additional  parameter  HRCON  for  each  material. 
HRCON  is  the  hourglass  control  factor  for  quadrilaterals  (0.05  to  0.2  is  recommended). 

Misoelianeous  Card  (716)  -  Add  one  additional  parameter  NTGRP  in  column  SO-^O. 
NTGRP  is  the  number  of  groups  into  which  the  finite  element  mesh  is  subdivided. 

Element  Data  Cards  (1016)  -  Add  one  additional  parameter  NETYPE  in  column  40-50. 
NETYPE  is  set  to  one  for  quadrilaterals,  otherwise  equal  to  sero  or  blank. 

Spedal  Shapes  Element  Data  Cards  (816)  -  Add  one  additional  parameter  NETYPE  (0  or  1) 
in  column  36-40. 

Element  Group  Array  for  Subcycling  (616)  -  The  parameters  are  IGRP,  KG  I,  KG2,  KGS 
and  KG4.  For  each  element  group  there  should  be  at  least  one  card.  This  group  of  cards  is 
followed  with  a  blank  card  to  flag  that  all  element  groups  are  completed.  This  group  of  cards 
is  inserted  right  after  all  element  cards  and  blanks.  IGRP  is  the  number  of  the  element 
groups  whereas  KGl,  KG2  and  KGS  are  indices  of  DO  loops  and  KG4  is  an  interval.  This 
can  best  be  explained  with  an  example.  Consider  the  elements  1  to  20  and  that  they 
represent  two  element  groups.  The  first  being  1,2  6  ,7  11,12  16,17  and  the  second  3,4,5 
8,0,10  13,14,15  18,10,20.  The  first  element  group  would  have  KG1=1  and  KG2=2  where 
1  and  2  are  the  first  and  last  entry  in  the  first  set.  KG 4 =5  is  the  interval  between  1  and  6 
and  KG3— 3  the  number  of  sets  less  one.  The  second  element  group  would  have  KG  1=3 
and  KG 2=5  and  an  interval  of  KG 4=5  (difference  of  3  and  8)  and  again  KG 4=3. 

Rigid  Interface  Card  (16) 

This  card  follows  the  regular  slideline  cards 

column  variable  name  deacription 

1-5  NRSL  number  of  rigid  interfaces 

to  be  read  in. 

Number  of  Primary  and  Secondary  Nodes  (216) 
column  variable  name  deacription 


1  -  5 


NMNO 


number  of  primary  nodes  for  this 
rigid  interface 


0-  10 


NSNO 


number  of  secondary  nodes  for 
this  rigid  interface 


Primary  Nodes  for  Ri^d  Interface  (1015) 


column  variable  name 


detcrtpHon 


IMNO(l)  primary  node  numbers 


6-  10 


IMNO(2) 


Secondary  Nodes  for  lUgid  Interface  (1015) 


column  variable  name 


description 


ISNO(  1)  secondary  node  numbers 


6-  10 


ISNO(2) 


These  last  three  sets  are  repeated  NRSL  number  of  times. 


» 


«.t « A 
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